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Abstract 

We present a double precision routine in Fortran for the precise and fast numerical evalu- 
ation of the two Master Integrals (Mis) of the equal mass two-loop sunrise graph for arbitrary 
momentum transfer in d — 2 and d — 4 dimensions. The routine implements the accelerated 
power series expansions obtained by solving the corresponding differential equations for the 
Mis at their singular points. With a maximum of 22 terms for the worst case expansion a 
relative precision of better than a part in 10 15 is achieved for arbitrary real values of the 
momentum transfer. 
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PROGRAM SUMMARY 



Title of program: sunem 
Version: 1.0 
Release: 1 
Catalogue identifier: 

Program obtainable from: http://www-ttp.physik.uni-karlsruhe.de/Progdata/ 

Computers: all 

Operating system: all 

Program language: F0RTRAN77 

Memory required to execute: Size: 1532k 

No. of bits in a word: up to 32 

No. of processors used: 1 

No. of bytes in distributed program: 34589 

Distribution format: ASCII 

Other programs called: none 

External files needed: none 

Keywords: multi-loop Feynman integrals, differential equations 

Nature of the physical problem: numerical evaluation of the two Master Integrals of the equal mass 
two-loop sunrise Feynman graph for arbitrary momentum transfer in d = 2 and d — 4 dimensions. 

Method of solution: accelerated power series expansions obtained by solving the differential equa- 
tions for the Mis at their singular points. With a maximum of 22 terms for the worse case expansion 
a relative precision of better than a part in 10 15 is achieved for arbitrary real values of the mo- 
mentum transfer. 

Restrictions on complexity of the problem: limited to real momentum transfer and equal internal 
masses. 

Typical running time: approximately 1 /j,s to evaluate the four Master integrals for a fixed momen- 
tum transfer value on a Pentium IV/3 GHz Linux PC. 



1 



LONG WRITE-UP 



1 Introduction. 

The continuous progress in the evaluation of radiative correections to various elementary particle 
processes, driven by the increasing number of current and future precise experimental results, is 
faced with the problem of the fast and precise numerical evaluation of those radiative corrections 
that cannot be expressed in terms of known analytic functions for which numerical routines are 
already available. 

In this paper we discuss a FORTRAN routine for the precise and fast numerical evaluation of the 
two Master Integrals (Mis) associated with the two-loop selfmass sunrise graph for equal masses, 
based on the results of pQ, where their analytic properties were established by studying the linear 
system of first order inhomogeneous differential equations satisfied by the Mis themselves. 

As a first obvious step, the linear system can be rewritten as a second order (inhomogeneous) 
linear differential equation for one of the Mis (from now on referred to as the main MI, or simply 
the MI), while the other MI can be expressed in terms of the first and its derivative. The accurate 
investigation of the second order equation shows that the immaginary part of the MI (which 
satisfies the associated homogeneous equation) is an ellyptic function of complicated argument (a 
fact known however since a long time), while the real part is likely to belong to a new family of 
analytic functions, which could not be traced back to the usual Nielsen's polylogarithms nor to 
their generalizations [2]. 

Nevertheless, the differential equation can be used to obtain a convenient integral represen- 
tation of the solution and, among other things, the explicit limiting behaviour of the MI under 
investigation at zero momentum transfer, at the pseudothreshold, at the threshold and at asymp- 
totic momenta. Given those limiting values as starting points, the differential equation can then 
be used to generate the coefficients of the relevant expansions up to any required order. 

In this paper we show that those expansions can be used for writing a (double precision) 
FORTRAN routine for the fast and precise numerical evaluation of both the Mis, within the d- 
continuous regularization scheme, in d = 2 and d = 4 dimensions. According to p^, it is indeed 
sufficient to work out the relevant expansions only for the main MI in d = 2 dimensions by direct 
iterative solution of the differential equation; the other MI in d = 2 dimensions, as well as both 
Mis in d = 4 dimensions, are then expressed in terms of the first MI in d = 2 dimensions and 
its derivatives (needless to say, differentiation becomes a simple algebraic operation when dealing 
with expansions). 

In the next 3 Sections we discuss the relevant expansions of the first MI in d = 2 dimensions 
for physical momentum around zero, at threshold and at infinity. In Section 6 we describe how to 
use the results of the previous Section for implementing the FORTRAN routine for the evaluation of 
the two Mis at d = 2 and d = 4 dimensions. In order to reach a relative precision of better than 

1 x 10 -15 (which is essentialy the limit of double precision FORTRAN) for any value of the squared 
momentun transfer s in the whole real axis, — oo < s < +oo, it is sufficient to keep a maximum of 
22 terms in the various expansions. 

2 Definition and notation. 

The two Mis associated with the (equal mass) 2-loop sunrise graph depicted in Fig. Q] are defined, 
in the d-continuous regularization scheme, as 




(1) 
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Figure 1: The 2-loop sunrise graph. 



where 



p is the external momentum, p 2 = z 
negative). One has, for arbitrary d 



V d k = 



1 d d k 



(2) 
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-s (when p is Euclidean or spacelike, z is positive and s is 



Si(d,z) = - 



-(d-3)+z 



dz 



S{d, z) 



At d = 2 S(d, z) and S\{d, z) are both finite, 

S(d,z) = 5(2, z) + 0{d - 2) , 
Sx(d,z) = S 1 (2,z) + 0(d-2) , 



(3) 



(4) 



and from now on we will simply consider 5(2, z), 5i(2, z). 5(2, z) satisfies the second order inho- 
mogenous equation 



dz 2 
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5(2, z) 



8z(z+ l)(z + 9) ' 
(5) 



which will be used to derive all the required expansions. 

Around d = 4 S(d, z) develops a double pole in [d — 4), so that it can be written as 

S(d, z) = j-J-^S^ (4, z) + (4, z) + S<°) (4, z) + 0(d - 4) . 



From pQ one has 



(d-4) 

5(- 2 )(4,z) 
5^(4, z) 



(6) 
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32 



(z + 18) , 



5W(4,z) = -L( z + l)( z + 9)(l + (z-3)^5(2,z) 



128 



(72 + 13z) 



(7) 



From Eq.©, recalling that 5i(2, z) is finite while Si(d, z) can be expanded around d = 4 as S(d, z) 
Eq.@, we find 



Si(2,z) 



1 + z- 



dz 



5(2, z) , 



(8) 
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S{- 2) (A,Z) 
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16 



-1 



dz 



S (°) (4, z)--(z+ 18) 



(9) 



We will work out explicitly, from the differential equations, the expansions of S(2, z) in the various 
regions; we will then use Ea.s (|8l7l9|) for obtaining the expansions of S\ (2, z), (4, z) and s[ ^ (4, z) 
in terms of the corresponding expansions of S(2, z) . 



3 The expansion around z = 0. 

The point z = (or s = 0) is a singular point of the differential equation Eq.ljSjl: the investigation 
carried out in pQ shows that the most general solution G{z) consists of a regular part and a singular 
part behaving like In z at z = 0, so that its expansion is 



oo oo 



G(z) = 2^ a n z n + In z ^ b n z n . (10) 

n=0 n=0 

Once the expansion Ea. p()|l is substituted in Eq.©, the differential equation takes the form of 
a system of linear algebraic equations for the two sets of coefficients a„, b n , and can be solved 
recursively expressing the coefficient of order n in terms of the coefficients of lower order. The 
procedure requires two initial conditions, i.e. the two coefficients ao and bo- 
As discussed in [J in the case of S(2, z) the two initial conditions are 

ao = ^?C1 2 = 0.146494... , 



12 V3, 

bo = 0, (11) 

where Q2 (f) is the Clausen function of argument 7r/3, while the vanishing of bo is due to the 
regularity of S*(2, z) at z = 0. 

Therefore, the expansion for S(2, z) reads simply 

00 

5(2,z) = ^a„z" , (12) 

n=0 

and on account of the explicit value of ao provided by Eq. (|llf) the differential equation then gives 
for the higher order coefficients 

V 7 ^, /7T\ 1 



at = -^Cl 2 (|) + - = -0.00716473..., 

K /7T\ 23 

^ Cl2 (3)- 864 -°- 000508185 - 



«2 
«3 



Vi^—Ch (-) + ^— = -0.0000414414... , (13) 
2916 V 3/ 7776 v ; 



and so on. The rational fractions (such as 1/24 in ai) are the expansion of the solution of the 
complete inhomogenous equation Eq.|0 with initial condition S*(2,0) = 0, whereas the irrational 
part with CI2 (^j, which is proportional to the regular solution of the associated homogenous 
equation, is fixed by the initial conditions. 

The nearest singular point of the differential equation Eq.© is at z = —1 (s = 1), so one would 
expect a convergence radius 1 for the expansion Eci. (|12fl : but the nearest singularity of S(2,z) 
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is at the physical threshold s = 9 or z = —9, so that 5(2, z) is analytic at the pseudothreshold 
z — — 1; that is achieved by suitable cancellations among the two (rational and irrational) terms 
which appear in the explicit expressions of the coefficients a n , as already shown also by the first 
coefficients Ea. (11811 . 

The next singular point of the differential equation is at z — —9 (or s — 9), which is also a non 
analyticity point for 5(2, z), as commented above; the expansion Eq. l|12f) has therefore convergence 
radius 9. 

The knowledge of the position of the nearest singularity allows to greatly speed up the conver- 
gence of the series, by introducing the "Bernoulli" variable (first proposed in [2] for the evaluation 
of Euler's dilogarithm, and repeatedly used in 4 for the evaluation of Harmonic Poly logarithms), 
which is given in this case by 

t = In (l + f). (14) 

One has 



z = 9(c £ -l)=9^+- , ( . 
replacing z in Eq. (|12|1 by its expansion in £, one obtains a new expansion in t which we write as 

oo 

S(2,z) = J2 a nt n ■ (15) 

n=0 

The explicit knowledge of N terms of the expansion Ea. l|12J) detrmines the first N terms of Ea. 1(15(1 . 
which for a same number of expansion terms and small t (see Section 6 for more details) approx- 
imates 5(2, z) much better than the original expansion in z. A straightforward calculation gives 
for the first coefficients 

a = a Q = j^Ch (J) = 0.146494... , 
ai = | - ^j-Ch (|) = -0.0644826... , 

a 2 = V3^C1 2 HT) - !!! = 0.00892173.. 

8 V 3 / 

a 3 = 



32 

-^^Cl 2 (|) + f = 0.000205110... , (16) 

and so on. Note the strong cancellation between the two terms in each coefficient of the expansion; 
the cancellation, already present in Ea. 1)13(1 as remarked above, is even bigger in Ea.l|T?))l. and 
grows quickly with the order n of the coefficient a n . Keeping the control of cancellations in the 
numerical evaluation of the coefficients would be in principle a delicate task, but that is no longer 
problem when using the arbitrary precision features of Mathematica 5 . 



4 The expansion around z = —9. 

As already remarked, the point z = —9 (s = 9) is the nearest singular point of 5(2, z). For z — * —9, 
the general solution consists of a regular part and a singular part behaving like ln(z + 9), so that 
5(2, z) admits an expansion of the general form 

oo oo 

5(2, z) = J2 a n (z + 9)" + ln(z + 9) ^ b n (z + 9)™ , (17) 

n=0 n=0 

where the coefficients, called again a n , b n as in Ea. (|10f> . are of course different from the coefficients 
of the expansion at z = 0. For z > —9, i.e. s < 9 or s below threshold, both terms of Eq. l|17|) are 
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real, while for z < —9, or s > 9, above threshold, ln(z + 9) develops an immaginary part, which 
with the usual s + ie prescription is —in, 



ln(> + 9) = ln|z + 9| -i n6(-z- 9) . (18) 
From we have the initial conditions 



°° = "IT 



°Q = -77T 



7rln(72)-5Cl 2 (|)] . (19) 



.3 

From the differential equation we found for the first coefficients 

at = ^ + b ^ + ^, 

h = b Q — , 

7 97 5 

° 2 " a °864 + bo 10368 + 6912 ' 

62 = 6o 864< ^ 

and so on for the higher order coefficients. As in the previous section, the terms independent of 
the initial conditions ao,6o (such as the fraction 1/192 in ax) are due to the inhomogenous term 
in the differential equation Eq.0. 

The nearest singular point for the two components of the expansion Ea. 1)17(1 is at z = —1, so 
that the radius of convergence of the two expansions is 8, i.e. the expansions converge in the range 
— 17 < z < — 1 or 1 < s < 17. That might look surprising, as 5(2, z) is analytic at z = — 1; but the 
two terms in the expansion Ea. l|17|l cannot combine to a single expansion converging in a wider 
range, given the presence of the \n(z + 9) in front of the second term. 

In this case the "Bernoulli" variable is therefore 

t = -ln(l-i±^, (21) 
and the expansion of S(2, z) in terms of t reads 

CO oo 

5(2, ») = £ a n t n + ln(z + 9) ^ f3 n t n . (22) 



n=0 n=0 



The explicit calculation gives 

a = a = 0.301695... , 
A, = b a = -0.113362... , 

aj = -a + -6 + 24 = 0.179817... , 

01 = ^b = -0.0755749... , 

5 26 11 
a 2 = —a + ^fco + ^ = 0.0449445... , 

f3 2 = -^b Q = -0.0209930... , (23) 
with the by now usual strong cancellations within the various terms of each coefficient. 
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5 The expansion at infinity. 

The differential equation has also z = oo as singular point; for z — ► oo the expected behaviour for 
5(2, z) is 



I ( °° I °° 1 °° 1 \ 

\n=0 n=0 n=0 / 



(24) 



For 2; > {i.e. s spacelike) the previous expression is real; in the timclikc region, we have to 
replace z by z — ie; that gives 

lnz = In \z\ — i7T0(— z) , 
ln 2 z = In 2 \z\ - 6{-z) (2i7iTn|z| +tt 2 ) , (25) 

from which one can obtain separately the real and immaginary parts of 5(2, z). 
From [J the initial conditions are 

ao = , 
b = 0, 

^ = y & ■ (26) 

It is to be noted that the presence of the term 3/16 In 2 z is anyhow required by the inhomogenous 
term in the differential equation, so that strictly speaking the initial conditions are simply ao = 
&o = 0. The differential equation then gives for the first coefficients 



a i 
h 
c\ 

(>2 
b 2 



3 

2 ' 
9 

~16 ' 
3 

^32 ' 
39 

~T ' 

45 



C2 - 16' (2?) 



and so on for the higher orders. 

As the nearest singularity is at z = —9, the "Bernoulli" variable is 



t = In I 1 + - ] , (28) 



z 



and we write the expansion of 5(2, z) in terms of t as 

/ oo oo oo \ 

5(2, z) = [J2 a nt n + \nzY / Pnt n + In 2 z £ 7n t" , (29) 



so that the first coefficients are 



a Q = p = 70 = 
ai = a — , 

/?! = = , 
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«2 

02 
72 
«3 
03 
73 



1 

48 ' 
1 

216 ' 
1 

54 ' 
1 

288 ' 
35 

776 ' 

5 

972 ' 
1 

2592 



(30) 



etc. 



6 The FORTRAN numerical routine. 

As a first step, we obtained the analytic expression of the coefficients of all the needed expansions 
for S(2,z) up to the required order (see below). From those we obtained the corresponding ex- 
pansions for the other three quantities, Si (2, z), S^(4, z) and S[ (4, z) by means of the formulae 
Eq.s l8l7l9ll . which involve only elementary algebra and differentiation. When dealing with power 
series expansions, differentiating is trivial, even if one term is lost at each differentiation; but that 
is easily compensated by allowing for a few extra terms in the initial expansion. 

We then evaluated the numerical value of all the coefficients so obtained. The analytic and 
numerical evaluations were carried out by Mathematica, in arbitrary precision arithmetics mode 
for keeping rounding errors under check. The numerical values of the coefficients were transferred 
to the FORTRAN code as constants with 18 digits; that exceeds of the course the double precision of 
FORTRAN but eliminates rounding problems in the least significant digit. 

In the range — oo < z < —17.45, or 17.45 < s < +oo (timelike region above threshold, all Mis 
complex), we evaluate the 4 functions by means of truncated expansions of the form 

N N N 

OLnt n + hi z ^ 0nt n + In 2 z £ ln t n , (31) 

n— n=0 n— 

where t is defined by Ea. (|28() and varies therefore in the range —0.725173... < t < 0, the complex 
value of In z is specified in Eg. (|25|) and the values of the coefficients a n ,0 n ,Jn depend of course 
on the Mis. We obtain the required double precision (relative error less than 1 x 10~ 15 ) for 
all the 4 functions with TV = 22. Those numerical checks, beyond the possibilities offered by 
double precision FORTRAN, were carried out by using Mathematica in arbitrary precision mode and 
including a growing number of terms in the expansion. 

In the range —17.45 < z < —5.15, or 5.15 < s < 17.45 (timelike region across the threshold 
s = 9) we use for the 4 functions truncated expansions of the form 

R N N 

— - + V« n i"+ln(z + 9) Y l 0nt n , (32) 
z + 9 ^-^ 

n=0 n=0 

where R — — \/3 7r/48 for S\{2,z) and vanishes in all the other cases, t, defined in Ea. (|21(l . varies 
in the range -0.720884... < t < 0.656333... and the value of ln(z + 9) is given in Ea.ljT3|). The 
required precision (better than 1 x 10~ 15 ) is obtained with N = 21. 
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In the range —5.15 < z < 11, or —11 < s < 5.15, we use truncated expansions of the form 

N 

^>„f\ (33) 

n=0 

where t, defined in Ea. (|14|l . varies in the range —0.849151 < t < 0.798508, and the required 
precision is obtained with N = 18. Note that the expansion in t is valid, and still quickly convergent, 
for corresponding values of z also in the region 9 < z < 11, which lies outside the region of 
convergence of the expansion in z from which the expansion in t was derived! 

Finally, in the range 11 < z < +oo, or — oo < s < —11 (spacclikc region), we use again the 
expansion Ea. (|31|) . with t varying in the range < t < 0.597837...; in this case lnz is real, and the 
aimed precision is obtained with N = 20. 

The numerical values of the 4 functions, evaluated according to the previous formulas, are 
returned by calling the FORTRAN routine 

SUBROUTINE SUNRISE (S , SR2 , SI2 , S1R2 , S1I2 , SR4 , SI4 , S1R4 , SI 14) , 
where all the variables are real DOUBLE PRECISION variables, S is the momentum transfer (S, equal 
to — z of previous sections, is positive when timelike, with physical threshold at 9), SR2,SI2 are 
the real and immaginary parts of S(2, z) as defined in Eq.QJ for d = 2 and z — — S , S1R2 , S1I2 the 
real and immaginary parts of Si (2, z) as defined also in Eq.JU, SR4,SI4 the real and immaginary 
parts of S^{A, z) as defined by Ea.s (jll6f ). and S1R4.S1I4 finally the real and immaginary parts of 
S , 1 °^(4, z). S can vary in the whole real axis, — oo < S < +oo; the immaginary parts, which appear 
when S is above threshold (S > 9), are evaluated according to the standard S + ie prescription. 

The calculation is carried out with relative precision better than 1 x 10 -15 (which is essentialy 
the double precision limit of FORTRAN), except of course when the corresponding quantities approach 
zero. 5(2, z) and S\{2,z) are positive definite and tend to zero for \z\ — » oo, while S^(4, z) and 
Si (4, z) have a zero around z — 12.5910... and z = —3.5599..., respectively. 

Internally, the routine SUNRISE calls, depending on the actual value of S, one of the three sub- 
routines SUNRISE0, SUNRISE9 and SUNRISE00. Those routines implement the various expansions 
as discussed in the previous sections. For ease of check and implementation, the routines contain 
more coefficients than actually used, each coefficient being an entry in a DATA statement written 
as a constant with 18 decimal digits (as already remarked above, that exceeds of course the double 
precision of FORTRAN but eliminates rounding problems). The arguments of the three routines 
consist of the same arguments of SUNRISE and an extra integer argument, specifying how many 
coefficients of the expansions are to be considered in the calculation. The three internal routines 
can in principle be called independently, for studying the contributions of the various terms of the 
expansions etc., as well as for additional checks. 

As a last remark, we succeeded in covering the whole real axis in z with just three different 
expansions and a maximum of 22 terms at worst; that is fully satisfactory, and from the practical 
point of view there is no need of further improvements. But let us observe that, in principle, the 
calculation could be further speeded up by using additional auxiliary expansions. For instance, 
one could consider another expansion point, say for instance at z = —18 (or s = 18); as that is a 
regular point in which the function is analytic, one can write as 

A + B{z + 18) + O {{z + 18) 2 ) 

the limiting behaviour of the function to be evaluated, parametrizing it in terms of two initial 
constants A, B, still to be determined. By using the differential equation, it is immediate to 
generate as many coefficients as desired of the expansion in powers of (z + 18). The two constants 
can then be obtained by matching numerically the new expansion with the expansion around 
z = —9 at, say, z = —13, and with the expansion at z = — oo at, say, z — 24. In so doing, one has 
to use every involved expansion on a smaller range, where the desired precision can be achieved 
with a smaller number of terms. 
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By the same token, one could consider the expansion around z = —9 Eg. l|17|l and keep the 
initial conditions ao,&0 7 whose exact value is given by Ea. (|19fl . as still unknown constants, use 
the differential equation to work out the desired coefficients of the expansion around z — — 9, and 
then obtain numerically ao,6o by imposing the smooth matching at, say, z — — 5 with the value 
obtained with the expansion around z = 0, supposedly known. One could then continue to recover 
in that manner the numerical value of the initial constants of the expansion at infinity, supposedly 
unknown as well. Note that in this procedure one can insert as many auxiliary expansions as 
considered useful; but one cannot skip the expansions at the singular points of the function, as 
those points cannot be included, by the very definition of singular points, in the convergence radius 
of any expansion around nearby points. The differential equation, again, gives the most general 
behaviour of the function at any point, to be used for obtaining the expansion around that point. 

7 Example 

The following simple program illustrates how to evaluate the real and imaginary parts of the 
Master integrals S(2,z), Si(2,z), S<°)(4, z) and sf ] (A,z) for a given value of the momentum 
transfer s — — z: 

program examplesunem 
implicit none 

real*8 S , SR2 , SI2 , S1R2 , S1I2 , SR4 , SI4 , S1R4 , S1I4 

write (6 , *)' Input squared momentum transfer S:' 
read(5,*) S 

call sunrise (S,SR2,SI2,S1R2,S1I2,SR4,SI4,S1R4,S1I4) 

write(6,*) 'First Master integral in d=2 dimensions' 
write (6, 200) S,SR2 
write(6,201) S,SI2 

write (6 ,*)' Second Master integral in d=2 dimensions' 
write (6, 210) S,S1R2 
write(6,211) S,S1I2 

write(6,*) 'First Master integral in d=4 dimensions: finite parts' 
write (6, 400) S,SR4 
write(6,401) S,SI4 

write (6 ,*)' Second Master integral in d=4 dimensions: finite parts' 
write(6,410) S.S1R4 
write(6,411) S,S1I4 



200 


format ( ' 


Re S(2, ' ,1PE22.15, ' 


) 


' ,1PE22. 


.15) 


201 


format ( ' 


Im S(2, ' ,1PE22.15, ' 


) 


' ,1PE22. 


.15) 


210 


format ( ' 


Re SI (2, ' ,1PE22.15, 


') 


' ,1PE22. 


.15) 


211 


format ( ' 


Im SI (2, ' ,1PE22.15, 


') 


' ,1PE22, 


.15) 


400 


format ( ' 


Re S~(0) (4, ' ,1PE22. 


15,') = 


' ,1PE22, 


.15) 


401 


format ( ' 


Im S~(0) (4, ' ,1PE22. 


15,') = 


' ,1PE22, 


.15) 


410 


format ( ' 


Re S1~(0) (4, ' ,1PE22 


.15, ')= 


' ,1PE22, 


.15) 


411 


format ( ' 


Im S1~(0) (4, ' ,1PE22 


.15, ')= 


' ,1PE22. 


.15) 




end 
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